/*
 heat-spectral lib. solves time-dependent PDEs using SEM-NI method.
 Copyright (C) 2009, Altieri Michele and Stefanoni Andrea
 Copyright (C) 2008, Proverbio Alessandro and Radice David
 
 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.
 
 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 GNU General Public License for more details.
 
 You should have received a copy of the GNU General Public License
 along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef FUNCTIONALS_HPP
#define FUNCTIONALS_HPP 1

#include "macro.hpp"

/**
    \file functionals.hpp
    \brief Useful functionals

    This file contains some useful linear and non-linear functional
    used to compute the numeric errors or to build the force term
    from a given function.

    This functionals are to be used with the evaluate functions in
    the solvers APIs
    \see double spectral::Problem::evaluate(std::valarray<double> const & u) const

*/

namespace spectral {

// ******************************************************************
// Standard Functionals
// ******************************************************************
// Some usefull functionals:
// 1. Standard linear functional
// 2. Time dependent linear functional
// 3. L2 Error functional (squared)
// 4. Time dependent L2 Error functional (squared)
// 5. H1 Error functional (squared)
// 6. Time dependent H1 Error functional (squared)
// 7. L2 Norm functional (squared)
// 8. H1 SemiNorm functional (squared)
// ******************************************************************

// ------------------------------------------------------------------
// Basic linear functional
// ------------------------------------------------------------------
/*! \brief Linear functional
 *
 *  If \e f is a functor of class \e FUNC, then
 *  the Functional \e F is such that \f$ F(v)=(f,v)_\Delta \f$.
 */
template<class MESH, class FUNC, bool TIME=false>
class Functional {
    public:
        /// The constructor
        Functional(MESH const & me, double const & t=0): me(me) {}

        /// Computes F(v)
        template<class V> inline
        double operator()(V const & v, int n) const {
            return me.prodL2(me.restriction(f,n),v,n);
        }
    private:
        MESH const & me;        /// mesh
        FUNC f;                 /// functor
};
// ------------------------------------------------------------------


// ------------------------------------------------------------------
// Time dependent linear functional
// ------------------------------------------------------------------
/*! \brief Time dependent linear functional
 *
 * If \e f(t) is a functor of class \e FUNC, then
 * the Functional \e F is such that \f$ F(v,t)=(f(t),v)_\Delta \f$.
 */
template<class MESH, class FUNC>
class Functional<MESH,FUNC,true>{
    private:
	struct ActualTime {
	    double t;
	    FUNC f;
	    ActualTime(double const & t) : t(t) {}
	    double operator()(double const & x, double const & y) const {
		return f(x,y,t);
	    }
	};
    public:
	Functional(MESH const & me, double const & t=0): me(me), f(t) {}

	template<class V> inline
	double operator()(V const & v, int n) const {
	    return me.prodL2(me.restriction(f,n),v,n);
	}
    private:
	MESH const & me;
	ActualTime f;
};

// ------------------------------------------------------------------
// L2 Error
// ------------------------------------------------------------------
/*! \brief Quadratic functional used to compute the \b squared L2 
 *         distance between a given functor and \e EXACTSOL
 */
template<class EXACTSOL,class MESH,bool TIME=false>
class L2Error {
    public:
        /// The constructor
        L2Error(MESH const * me): me(me) {}

        /// Computes F(v)
        template<class T>
        double operator()(T const & unum, int n) const {
            return
                SPECTRAL_RE_Q_FORM(unum,uex);
        }
    private:
        MESH const * me;        ///< mesh
        EXACTSOL uex;           ///< exact solution
};

// ------------------------------------------------------------------
// Time dependent L2 Error
// ------------------------------------------------------------------
/*! \brief Quadratic functional used to compute the \b squared
 * 	   L2 distance between a given functor and
 * 	   \e EXACTSOL at a fixed time \b t.
 *  \param EXACTSOL the exact solution
 *  \param t the actual time passed to the constructor
 *  In particular, if \f$ u \f$ is the exact solution of the problem
 *  and \f$ u_n \f$ is the numerical one, this function computes the
 *  error norm
 *  \f[
 *  	\|u(t)-u_n(t)\|^2_{L^2(\Omega)}
 *  \f]
 */
template<class EXACTSOL,class MESH>
class L2Error<EXACTSOL,MESH,true> {
    private:
	template<class FUN>
	struct ActualTime {
	    double t;
	    FUN u;
	    ActualTime(double const & t) : t(t) {}
	    double operator()(double const & x, double const & y) const {
		return u(x,y,t);
	    }
	};

    public:
	L2Error(MESH const * me, double const & t): me(me), uex(t) {}

	template<class T>
	double operator()(T const & unum, int n) const {
	    return
		SPECTRAL_RE_Q_FORM(unum,uex);
	}

    private:
	MESH const * me;
	ActualTime<EXACTSOL> uex;
};
// ------------------------------------------------------------------

// ------------------------------------------------------------------
// H1 Error
// ------------------------------------------------------------------
/*! \brief Quadratic functional used to compute the \b squared H1 
 *         distance between a given functor and \e U
 *  \param U the exact solution
 *  \param UX the partial derivatives with respect to x of \e U
 *  \param UY the partial derivatives with respect to y of \e U
 */
template<class U,class UX,class UY,class MESH,bool TIME=false>
class H1Error {
    public:
        /// The constructor
        H1Error(MESH const * me): me(me) {}

        /// Computes F(unum)
        template<class T>
        double operator()(T const & unum,int n) const {
            return
                + SPECTRAL_RE_Q_FORM(unum,u)
                + SPECTRAL_OP_RE_Q_FORM(unum,ux,Dx)
                + SPECTRAL_OP_RE_Q_FORM(unum,uy,Dy);
        }
    private:
        MESH const * me;  ///< mesh
        U u;              ///< exact solution
        UX ux;         /**< \brief partial derivatives of the exact solution
                                   with respect to x */
        UY uy;         /**< \brief partial derivatives of the exact solution
                                   with respect to y */
};

// ------------------------------------------------------------------
// Time dependent H1 Error
// ------------------------------------------------------------------
/*! \brief Quadratic functional used to compute the \b squared
 * 	   H1 distance between a given functor and
 * 	   \e U at a fixed time.
 *  \param U the exact solution
 *  \param UX the partial derivatives with respect to x of \e U
 *  \param UY the partial derivatives with respect to y of \e U
 *  \param t the actual time passed to the constructor
 *  In particular, if \f$ u \f$ is the exact solution of the problem
 *  and \f$ u_n \f$ is the numerical one, this function computes the
 *  error norm
 *  \f[
 *  	\|u(t)-u_n(t)\|^2_{H^1(\Omega)}
 *  \f]
 */
template<class U,class UX,class UY,class MESH>
class H1Error<U,UX,UY,MESH,true> {
    private:
	template<class FUN>
	struct ActualTime {
	    double t;
	    FUN u;
	    ActualTime(double const & t) : t(t) {}
	    double operator()(double const & x, double const & y) const {
		return u(x,y,t);
	    }
	};

    public:
	H1Error(MESH const * me, double const & t): me(me), u(t), ux(t), uy(t) {}

	template<class T>
	double operator()(T const & unum,int n) const {
	    return
		+ SPECTRAL_RE_Q_FORM(unum,u)
		+ SPECTRAL_OP_RE_Q_FORM(unum,ux,Dx)
		+ SPECTRAL_OP_RE_Q_FORM(unum,uy,Dy);
	}
    private:
	MESH const * me;
	ActualTime<U> u;
	ActualTime<UX> ux;
	ActualTime<UY> uy;
};
// ------------------------------------------------------------------

// ------------------------------------------------------------------
// L2 Norm
// ------------------------------------------------------------------
/// Quadratic functional used to compute the \b squared L2 norm of a function
template<class MESH>
class L2Norm {
    public:
        /// Constructor
        L2Norm(MESH const * me): me(me) {}

        /// Computes F(v)
        template<class T>
        double operator()(T const & unum, int n) const {
            return
                me->prodL2(me->restriction(unum,n),
                           me->restriction(unum,n));
        }
    private:
        MESH const * me;        ///< mesh
};
// ------------------------------------------------------------------

// ------------------------------------------------------------------
// H1 SemiNorm
// ------------------------------------------------------------------
/*! \brief Quadratic functional used to compute the \b squared H1 
           semi-norm of a function
*/
template<class MESH>
class H1SemiNorm {
    public:
        /// Constructor
        H1SemiNorm(MESH const * me): me(me) {}

        /// Computes F(v)
        template<class T>
        double operator()(T const & unum,int n) const {
            return
                + me->prodL2(me->Dx(me->restriction(unum,n),n),
                             me->Dx(me->restriction(unum,n),n))

                + me->prodL2(me->Dy(me->restriction(unum,n),n),
                             me->Dy(me->restriction(unum,n),n));
        }
    private:
        MESH const * me;        /// mesh
};
// ------------------------------------------------------------------

// ------------------------------------------------------------------
// H1 Norm
// ------------------------------------------------------------------
/// Quadratic functional used to compute the \b squared H1 norm of a function
template<class MESH>
class H1Norm {
    public:
        /// Constructor
        H1Norm(MESH const * me): me(me) {}

        /// Computes F(v)
        template<class T>
        double operator()(T const & unum,int n) const {
            return
                + me->prodL2(me->Dx(me->restriction(unum,n),n),
                             me->Dx(me->restriction(unum,n),n))

                + me->prodL2(me->Dy(me->restriction(unum,n),n),
                             me->Dy(me->restriction(unum,n),n))

                + me->prodL2(me->restriction(unum,n),
                             me->restriction(unum,n));
        }
    private:
        MESH const * me;        /// mesh
};
// ------------------------------------------------------------------

} //namespace

#endif
